Coral health statuses were monitored for disease in the CCMI nursery Dataset is located in Dryad: Brown, Anya et al. (2022), CCMI nursery coral disease 2019, Dryad, Dataset, https://doi.org/10.25338/B8F643

library(tidyr)
library(dplyr)
library(reshape2)
library(plyr)
library(ggplot2)
library(Rmisc)
library(stringr)
library(lme4)
library(lmerTest)
library(cowplot)
library(patchwork)
library(viridis)
library(glmmTMB)
library(DHARMa)
#dataset in Dryad

Data set up

dis_1 <-read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Frame Surveys/Data/FrameSurvey_Dec2019_wideformat_edit2.csv")

setwd("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Frame Surveys/Data/")
#if goes from dead to missing, change back to dead 
library(data.table)

Make sure data in columns are consistently named

dis <- reshape2::melt(dis_1, id= c("Frame", "Column","Row","Genotype")) #value.name ="Coral_stat")
colnames(dis)[5] <- "Date"
colnames(dis)[6] <- "Coral_stat"

#remove X from the dates
dis$Date <-gsub("X","",dis$Date)

#create unique IDs for each coral based on frame number and location on a frame
dis$ID <- paste(dis$Frame, dis$Column, dis$Row)

#remove Xs from data
#dis$Date <- str_remove(dis$Date,"X")
#remove leading or trailing white spaces
dis$Coral_stat <- trimws(dis$Coral_stat)

#Replace P with H
p <- which(dis$Coral_stat == "P")
dis[p,"Coral_stat"] <- "H"


#Replace XS with X
s <- which(dis$Coral_stat == "XS")
dis[s,"Coral_stat"] <- "X"

#Replace HR with H (Healty recovred should be just health)
h <- which(dis$Coral_stat == "HR")
dis[h,"Coral_stat"] <- "H"

#replace "0" with O for missing
o <- which(dis$Coral_stat == "0")
dis[o,"Coral_stat"] <- "O"

#replace OD with Dead
d <- which(dis$Coral_stat == "OD")
dis[d,"Coral_stat"] <- "D"

#replace R for Recovered with Healthy
e <- which(dis$Coral_stat == "R")
dis[e,"Coral_stat"] <- "H"

#replace diseased and recovered (XR) with X (disease)
f <- which(dis$Coral_stat == "XR")
dis[f,"Coral_stat"] <- "X"

#replace recovered with disease with just disease
g <- which(dis$Coral_stat == "RX")
dis[g,"Coral_stat"] <- "X"

#remove data from May 16 - not trustworthy
dis <- dis[which(dis$Date != "2019.05.16"),]
#remove missing corals (never present) and where there is no data
dis <- dis[which(dis$Coral_stat != "O" & dis$Coral_stat != "NoData"),]

Genotype counts

counts_dis_g <- plyr::ddply(na.omit(dis), c("Date", "Genotype","Frame","Coral_stat"), summarize,
      count = length(Coral_stat))

Including zeroes for other health categories

#put all data longwise 
counts_dis_g <- as.data.table(counts_dis_g)
dim(counts_dis_g)
## [1] 2668    5
counts_cast <- dcast.data.table(counts_dis_g, Date+Genotype+Frame~Coral_stat, value.var = "count")

counts_cast <- counts_cast %>% mutate_all(~replace(., is.na(.), 0))
dim(counts_cast)
## [1] 1433    6

Melt data

counts_cast_melt <- melt(data = counts_cast, id.vars = c("Date","Frame","Genotype"))
colnames(counts_cast_melt)[4] <- "Coral_stat"
colnames(counts_cast_melt)[5] <- "count"

Labeling single and mixed genotype frames

library(data.table)
frame <- data.table(unique(counts_cast_melt$Frame))
colnames(frame)<- "Frame"
frame$sing_mix <- ifelse(frame$Frame > 17, "mixed","single")
s <- which(frame$Frame == "1207")
frame[s,"sing_mix"] <- "single"
w <- which(frame$Frame == "1279")
frame[w,"sing_mix"] <- "single"

Counts across all of the frames and genotypes

#count number of corals in a health category - this sums the number of corals over all health categories on a frame to get the total number of corals on a frame

counts_fr <- ddply(counts_cast_melt, c("Date", "Frame"), summarize,       countf = sum(count))

#for each frame, wheter it's a single or mixed genotype
singmixf <- merge(frame, counts_fr)

#total number of corals on single or mixed frames at each date
total.per <- ddply(singmixf, c("Date","sing_mix"), summarize, 
      total = sum(countf))

#total number of frames 
total.f <- ddply(singmixf, c("Date","sing_mix"), summarize, 
      total = length(Frame))
#12 mixed frames
#18 single frames 

#432 single corals
#218 mixed

Calculation of proportion of health statuses on a frame including zeroes

counts_all_singmix <- merge(counts_cast_melt, singmixf, by = c("Frame", "Date"))
dim(counts_all_singmix)
## [1] 4299    7
counts_all_of_frame <- ddply(counts_all_singmix, c("Date", "Frame","sing_mix","countf","Coral_stat"), summarize,       
                  count_all = sum(count))
dim(counts_all_of_frame)
## [1] 1872    6
#average over all frames
sum_all_counts <- Rmisc::summarySE(counts_all_of_frame, measurevar = "count_all", groupvars = c("Date","Coral_stat","sing_mix"))


counts_all_of_frame$prop <- counts_all_of_frame$count_all/counts_all_of_frame$countf

#average of proportion in each health category per frame 
sum_all <- Rmisc::summarySE(counts_all_of_frame, measurevar = "prop", groupvars = c("Date","Coral_stat","sing_mix"))

Figure 2: Plot of Disease for all corals across the nursery

sum_all2 <- subset(sum_all, Date != "2019.02.01")
sum_all2$Date2 <- gsub("[.]", "/", sum_all2$Date)
sum_all2$Date3 <- as.Date(sum_all2$Date2)

#### Used in Figure 2 below ####
sing_mix_all  <- ggplot(sum_all2[which(sum_all2$Coral_stat=="X"),], aes(x = Date3, y = prop, group =sing_mix)) + geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5)) + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom") + xlab("Date")

Disease-only data

disease <- counts_all_of_frame[which(counts_all_of_frame$Coral_stat == "X"),]
dim(disease)
## [1] 624   7
dead <- counts_all_of_frame[which(counts_all_of_frame$Coral_stat == "D"),]


total.disease.counts <- ddply(disease, c("Date","sing_mix"), summarize, 
sum.disease = sum(count_all), sum.all = sum(countf))

total.dead.counts <- ddply(dead, c("Date","sing_mix"), summarize, sum.disease = sum(count_all), sum.all = sum(countf))


#create factors and turn dates into date data type
disease$Frame <- as.factor(disease$Frame)
disease$Date2 <- gsub("[.]", "/", disease$Date)
disease$Date3 <- as.Date(disease$Date2)
#disease$Date2 <- as.numeric(disease$Date)

#includind a polynomial
x2 <-poly(disease$Date3,2)

Model 1 analysis (all data)

disease <- counts_all_of_frame[which(counts_all_of_frame$Coral_stat == "X"),]
disease$Frame <- as.factor(disease$Frame)
disease$Date2 <- gsub("[.]", "/", disease$Date)
disease$Date3 <- as.Date(disease$Date2)

#include density
Frame_area <- read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Frame Surveys/Data/Frame_area.csv")

#dim(disease)
#dim(disease_dens)
disease_dens <- merge(disease, Frame_area, by = "Frame")
disease_dens$density <- disease_dens$countf/disease_dens$Area
disease_dens$density_rescale <- reshape::rescaler(to = c(0,1), x = disease_dens$density)

Model 1: Binomial analysis

x2 <-poly(disease_dens$Date3,2)

library(glmmTMB)
glmmtb_alldat <- glmmTMB(cbind(count_all, countf-count_all) ~ sing_mix*x2+density + (1|Frame), family = binomial, data= disease_dens)

Check Model 1 residuals

library(DHARMa, quietly = TRUE)
glmtb_resid1 <- simulateResiduals(glmmtb_alldat)
testOutliers(glmtb_resid1, type = "bootstrap")

## 
##  DHARMa bootstrapped outlier test
## 
## data:  glmtb_resid1
## outliers at both margin(s) = 0, observations = 624, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01686699
## sample estimates:
## outlier frequency (expected: 0.00259615384615385 ) 
##                                                  0
plot(glmtb_resid1)

Summary of Model 1 as a table

summary(glmmtb_alldat)
##  Family: binomial  ( logit )
## Formula:          
## cbind(count_all, countf - count_all) ~ sing_mix * x2 + density +  
##     (1 | Frame)
## Data: disease_dens
## 
##      AIC      BIC   logLik deviance df.resid 
##   2997.4   3032.9  -1490.7   2981.4      616 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Frame  (Intercept) 2.333    1.527   
## Number of obs: 624, groups:  Frame, 30
## 
## Conditional model:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -1.6730     0.5349  -3.128  0.00176 ** 
## sing_mixsingle       1.0493     0.6604   1.589  0.11209    
## x21                 11.3317     2.1391   5.297 1.17e-07 ***
## x22                -23.3456     2.4903  -9.375  < 2e-16 ***
## density             -0.1371     0.1097  -1.250  0.21140    
## sing_mixsingle:x21  12.1196     2.6992   4.490 7.12e-06 ***
## sing_mixsingle:x22  -8.5804     3.1386  -2.734  0.00626 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(car::Anova(glmmtb_alldat))
Chisq Df Pr(>Chisq)
sing_mix 2.658957 1 0.1029682
x2 358.706944 2 0.0000000
density 1.561784 1 0.2114042
sing_mix:x2 23.471265 2 0.0000080

Density information

disease_dens %>%
   dplyr::group_by(Date,sing_mix)%>%
  dplyr::summarize(total = sum(count_all))
## `summarise()` has grouped output by 'Date'. You can override using the
## `.groups` argument.

Figure 1 data summary

sum_all2 <- Rmisc::summarySE(counts_all_of_frame, measurevar = "prop", groupvars = c("Date","Coral_stat"))

sum_all2$Date2 <- gsub("[.]", "/", sum_all2$Date)
sum_all2$Date3 <- as.Date(sum_all2$Date2)

Figure 1 Plot

#### Figure 1 ####

ggplot(sum_all2[which(sum_all2$Date!="2019.02.01"),], aes(x = Date3, y = prop, group = Coral_stat)) + geom_point(aes(fill = Coral_stat),color = "black", pch = 21, size = 3) + theme_bw() + geom_errorbar(aes(ymax = prop+se, ymin = prop-se)) + ylab("Proportion of health status over time")  + theme(axis.text.x = element_text(angle = 90), panel.grid.minor = element_blank(), panel.grid.major = element_blank()) + scale_fill_manual(values = c("black","orange", "gray"), name = "Coral Health", labels = c("Dead","Healthy", "Disease")) + xlab("Date")

Figure S2: Health statuses across all frames set up

#### Supplement - all Frames
sum_all3 <- Rmisc::summarySE(counts_all_of_frame, measurevar = "prop", groupvars = c("Date","Coral_stat", "Frame", "sing_mix"),na.rm = T)
## Warning in qt(conf.interval/2 + 0.5, datac$N - 1): NaNs produced
sum_all3$Date2 <- gsub("[.]", "/", sum_all3$Date)
sum_all3$Date3 <- as.Date(sum_all3$Date2)

sum_all3$Frame <- factor(sum_all3$Frame, levels = c("1235","1231","14","9","1237","7","5","1221","1219","12","11","1251","1201","1217","13","10","1209","6","4","3","8","1279","1207","16","15","1107","1","2","1295","1296"))

Figure S2: Health statuses across all frames

ggplot(sum_all3[which(sum_all3$Date!="2019.02.01"),], aes(x = Date3, y = prop, group = Coral_stat)) + theme_bw() + geom_errorbar(aes(ymax = prop+se, ymin = prop-se)) + ylab("Proportion of health status over time")  + theme(text = element_text(size = 20), axis.text.x = element_text(size = 16, angle = 90), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background = element_rect(color = "black", fill ="white")) + scale_color_manual(values = c("black","orange", "gray"), name = "Coral Health",labels = c("Dead","Healthy","Disease")) + xlab("Date") + facet_wrap(~Frame+sing_mix) + geom_line(aes(color = Coral_stat), size = 1)+ geom_point(aes(color = Coral_stat), size = 1)

Including genotype information calculates the proportion of disease within a genotype on a frame

counts_cast_melt <- melt(data = counts_cast, id.vars = c("Date","Frame","Genotype"))

colnames(counts_cast_melt)[4] <- "Coral_stat"
colnames(counts_cast_melt)[5] <- "count"

counts_cast_melt <- subset(counts_cast_melt, Coral_stat == "X" | Coral_stat == "H" | Coral_stat == "D")
dim(counts_cast_melt)
## [1] 4299    5
#Blue
B <- counts_cast_melt[which(counts_cast_melt$Genotype == "B"),]
counts_fr <- ddply(B, c("Date", "Frame"), summarize,      
                   countf = sum(count))
Bf <- merge(B, counts_fr, by = c("Date","Frame"))

Bf$prop <- Bf$count/Bf$countf            

genotypefunc <- function(Gen){
   Gent <- counts_cast_melt[which(counts_cast_melt$Genotype == Gen),]
   counts_fr <- ddply(Gent, c("Date", "Frame"), summarize,      
                   countf = sum(count))
Gentf <- merge(Gent, counts_fr, by = c("Date","Frame"))
Gentf$prop <- Gentf$count/Gentf$countf            
Gentf <- as.data.frame(Gentf)
}

Rf <- genotypefunc(Gen = "R")
dim(Rf)
## [1] 615   7
Kf <- genotypefunc(Gen = "K")
dim(Kf)
## [1] 501   7
Gf <- genotypefunc(Gen = "G")
dim(Gf)
## [1] 678   7
Yf <- genotypefunc(Gen = "Y")
dim(Yf)
## [1] 678   7
#combine together
GenRYGBK <- rbind(Rf, Kf, Bf, Gf, Yf)
dim(GenRYGBK)
## [1] 2970    7

Summarizing the data with genotypes

#merge with single mix

colnames(singmixf)[4] <- "count_all"
GenRYGBK_sm <- merge(singmixf, GenRYGBK, by = c("Date","Frame"))
#dimension check
dim(singmixf)
## [1] 624   4
dim(GenRYGBK)
## [1] 2970    7
dim(GenRYGBK_sm) #summarize - just B, K, R, Y, G
## [1] 2970    9
sumgen1 <- Rmisc::summarySE(GenRYGBK_sm, measurevar = c("prop"), groupvars = c("Date","Genotype","sing_mix", "Coral_stat"))

sum_onlygrybk <- Rmisc::summarySE(GenRYGBK_sm, measurevar = c("prop"), groupvars = c("Date","sing_mix", "Coral_stat"))

Model 2: Including genotypes in understanding disease prevalence

disease2 <- GenRYGBK_sm[which(GenRYGBK_sm$Coral_stat == "X"),]
disease2$Frame <- as.factor(disease2$Frame)
disease2$sing_mix <- as.factor(disease2$sing_mix)
disease2$Genotype <- as.factor(disease2$Genotype)

disease2$Date2 <- gsub("[.]", "/", disease2$Date)
disease2$Date3 <- as.Date(disease2$Date2)
#between genotype and single mixed
disease2_dens <- merge(Frame_area, disease2, by = "Frame")
dim(disease2)
## [1] 990  11
dim(disease2_dens)
## [1] 990  12
disease2_dens$density <- disease2_dens$countf/disease2_dens$Area
library(stats)
x2 <-poly(disease2_dens$Date3,2)
library("glmmTMB")

glmmtb_all <- 
glmmTMB(cbind(count, countf-count) ~ sing_mix*x2*Genotype + density + (1|Frame), family = binomial, data= disease2_dens)
summary(glmmtb_all)
##  Family: binomial  ( logit )
## Formula:          
## cbind(count, countf - count) ~ sing_mix * x2 * Genotype + density +  
##     (1 | Frame)
## Data: disease2_dens
## 
##      AIC      BIC   logLik deviance df.resid 
##   2622.6   2779.4  -1279.3   2558.6      958 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Frame  (Intercept) 1.611    1.269   
## Number of obs: 990, groups:  Frame, 26
## 
## Conditional model:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -3.01091    0.55771  -5.399 6.71e-08 ***
## sing_mixsingle                 2.76577    1.14032   2.425 0.015290 *  
## x21                           70.04911   18.12525   3.865 0.000111 ***
## x22                          -61.49875   18.21510  -3.376 0.000735 ***
## GenotypeG                      0.09659    0.31877   0.303 0.761891    
## GenotypeK                      0.87877    0.37570   2.339 0.019335 *  
## GenotypeR                      0.62232    0.31792   1.957 0.050293 .  
## GenotypeY                      0.72232    0.32005   2.257 0.024013 *  
## density                       -0.06856    0.12141  -0.565 0.572262    
## sing_mixsingle:x21           -53.44289   18.56310  -2.879 0.003990 ** 
## sing_mixsingle:x22            27.44410   18.82341   1.458 0.144847    
## sing_mixsingle:GenotypeG      -2.37924    1.01933  -2.334 0.019589 *  
## sing_mixsingle:GenotypeK      -1.00748    1.04460  -0.964 0.334813    
## sing_mixsingle:GenotypeR      -1.16649    1.09148  -1.069 0.285193    
## sing_mixsingle:GenotypeY      -1.01770    1.11115  -0.916 0.359721    
## x21:GenotypeG                -42.77861   19.46394  -2.198 0.027961 *  
## x22:GenotypeG                 54.64573   19.68237   2.776 0.005497 ** 
## x21:GenotypeK                 39.67457   24.72132   1.605 0.108522    
## x22:GenotypeK                -35.15077   25.40060  -1.384 0.166403    
## x21:GenotypeR                -79.63663   18.81900  -4.232 2.32e-05 ***
## x22:GenotypeR                 33.73691   18.98273   1.777 0.075528 .  
## x21:GenotypeY                -54.68006   19.79114  -2.763 0.005730 ** 
## x22:GenotypeY                 20.12562   20.57758   0.978 0.328056    
## sing_mixsingle:x21:GenotypeG  42.88986   20.70042   2.072 0.038272 *  
## sing_mixsingle:x22:GenotypeG -48.01809   21.37918  -2.246 0.024703 *  
## sing_mixsingle:x21:GenotypeK  19.36679   25.44816   0.761 0.446640    
## sing_mixsingle:x22:GenotypeK   6.29700   26.26124   0.240 0.810498    
## sing_mixsingle:x21:GenotypeR  63.39830   19.74293   3.211 0.001322 ** 
## sing_mixsingle:x22:GenotypeR -50.79416   20.61550  -2.464 0.013744 *  
## sing_mixsingle:x21:GenotypeY  63.16198   20.82326   3.033 0.002419 ** 
## sing_mixsingle:x22:GenotypeY -29.32972   21.98987  -1.334 0.182275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#anova(glm3a)
#car::Anova(glmmtb_all, type = "III")

Model 2 summary table

knitr::kable(car::Anova(glmmtb_all, type = "III"))
Chisq Df Pr(>Chisq)
(Intercept) 29.1458553 1 0.0000001
sing_mix 5.8827825 1 0.0152896
x2 15.0168983 2 0.0005484
Genotype 16.2412186 4 0.0027120
density 0.3189123 1 0.5722620
sing_mix:x2 15.1217178 2 0.0005204
sing_mix:Genotype 5.8737111 4 0.2087802
x2:Genotype 109.9652765 8 0.0000000
sing_mix:x2:Genotype 39.1258261 8 0.0000047 
library(DHARMa, quietly = TRUE)
glmtb_resid <- simulateResiduals(glmmtb_all)
testOutliers(glmtb_resid, type = "bootstrap")

## 
##  DHARMa bootstrapped outlier test
## 
## data:  glmtb_resid
## outliers at both margin(s) = 0, observations = 990, p-value = 0.72
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.008611111
## sample estimates:
## outlier frequency (expected: 0.00237373737373737 ) 
##                                                  0
plot(glmtb_resid)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

#no significant outliers

Figure 2

disease2 <- disease2[which(disease2$Date!= "2019.02.01"),]
sum_gen_all_sing_mix2 <- Rmisc::summarySE(disease2, measurevar = c("prop"), groupvars = c("Date3","sing_mix"), na.rm = T)

#summary for Figure 2
sum_gybrk_sing_mix <- Rmisc::summarySE(na.omit(disease2), measurevar = c("prop"), groupvars = c("Date","Date3","sing_mix","Genotype"))

#summary for SFig 3
sumall_gybrk_sing_mix <- Rmisc::summarySE(na.omit(GenRYGBK_sm), measurevar = c("prop"), groupvars = c("Date","sing_mix","Genotype","Coral_stat"))

Supplement Figure S3: Proportion of health status by genotypes on both single and mixed frames

sumall_gybrk_sing_mix$Date2 <- gsub("[.]", "/", sumall_gybrk_sing_mix$Date)
sumall_gybrk_sing_mix$Date3 <- as.Date(sumall_gybrk_sing_mix$Date2)
                          
sumall_gybrk_sing_mix$Coral_stat <- factor(sumall_gybrk_sing_mix$Coral_stat, levels = c("H","X","D"), labels = c("Healthy","Disease","Dead"))


ggplot(sumall_gybrk_sing_mix[which(sumall_gybrk_sing_mix$Date != "2019.02.01"),], aes(x = Date3, y = prop, group =sing_mix))  + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Proportion of Corals") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity")+ theme(text = element_text(size = 16), axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(vars(Genotype), vars(Coral_stat))

Individual genotype plots for disease prevalence

plot.Y <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="Y"),], aes(x = Date3, y = prop, group =sing_mix))  + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype) + lims(y = c(0,0.85))

plot.G <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="G"),], aes(x = Date3, y = prop, group =sing_mix))  + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))

plot.R <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="R"),], aes(x = Date3, y = prop, group =sing_mix))  + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))

plot.B <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="B"),], aes(x = Date3, y = prop, group =sing_mix))  + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))

plot.K <- ggplot(sum_gybrk_sing_mix[which(sum_gybrk_sing_mix$Genotype=="K"),], aes(x = Date3, y = prop, group =sing_mix))  + geom_errorbar(aes(ymax = prop+se, ymin = prop-se), position = position_dodge(width = 0.5), color = "gray") + ylab("Disease Prevalance") + theme_bw() + scale_fill_manual(values = c("darkorchid4","gold"), labels = c("mixed", "single"), name = "Diversity") + theme(axis.text.x = element_text(angle = 90),panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "bottom")+ geom_point(aes(fill = sing_mix), pch = 21, color = "black",size = 3, alpha= 0.75, position = position_dodge(width = 0.5)) + xlab("Date") + facet_grid(~Genotype)+ lims(y = c(0,0.85))

Figure 2

library(patchwork)

(sing_mix_all + plot_layout(widths = 25))/(plot.G + plot.Y + plot.R + plot.B + plot.K
+ plot_layout(nrow = 1, widths = 25))/plot_layout(guides = "collect") & theme(legend.position = "bottom") & plot_annotation(tag_levels = 'a')

Supplemental Figure 1: Temperature

tempcsv <- read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Microbial Sampling_wholeframes/Temperature/tempcsv.csv", header=FALSE, comment.char="#")

temp <- tempcsv[,2:3]
colnames(temp) <- c("Date_time","temp")
temp <- separate(temp, "Date_time",c("Date","time")," ")
temp$Date <- as.Date(temp$Date, format = "%m/%d/%y")
temp2 <- subset(temp, Date < "2019-08-19")
temp.sum <- summarySE(temp2, measurevar = "temp", groupvars = c("Date"))

martha_r <- read.csv("~/Dropbox/POST DOC/PostDoc Project/Spring 2019/Research/Research/Disease_deep nursery/Microbial Sampling_wholeframes/Temperature/martha_r.csv")

temp3 <-martha_r
colnames(temp3) <- c("Date_time","temp")
temp4 <- separate(temp3, "Date_time",c("Date","time")," ")
str(temp4)
## 'data.frame':    11179 obs. of  3 variables:
##  $ Date: chr  "2/11/19" "2/11/19" "2/11/19" "2/11/19" ...
##  $ time: chr  "15:30" "16:00" "16:30" "17:00" ...
##  $ temp: num  80.7 80.7 80.6 80.6 80.6 ...
temp4$Date <- as.Date(temp4$Date, format = "%m/%d/%y")

temp5 <- subset(temp4, Date > "2019-05-19" & Date < "2019-08-10")
temp.sum2 <- summarySE(temp5, measurevar = "temp", groupvars = c("Date"))

dim(temp.sum2)
## [1] 82  6
temp.sum2$Loc <- "Martha"
temp.sum$Loc <- "Deep"

temp.both <- data.table(rbind(temp.sum, temp.sum2))
temp.both$tempC <- (temp.both$temp - 32)* (5/9)
templab <- "Temperature (°C)"

templot <- ggplot(temp.both, aes(x = Date, y = tempC)) + geom_point(size = 3, fill = "black", color = "black", pch = 21) + scale_x_date(date_labels = "%m-%d",date_breaks = "10 day" ) + ylab(templab) + theme(text = element_text(size = 16)) + theme_bw() + geom_line()
library(cowplot)
sum_all2_a <- select(sum_all2, Coral_stat, N, prop, se, Date3)
colnames(sum_all2_a)[5] <- "Date"
temp.dis <- merge(sum_all2_a, temp.both, by = "Date")

tempdeadp <- ggplot(temp.dis[which(temp.dis$Coral_stat == "D"),], aes(x = tempC,  y = prop)) + geom_point(size = 3) + theme_bw() + xlab(templab) + ylab("Proportion of Dead Corals") + geom_smooth(se = F, col = "black")

tempdisp <-ggplot(temp.dis[which(temp.dis$Coral_stat == "X"),], aes(x = tempC,  y = prop)) + geom_point(size = 3) + theme_bw() + xlab(templab) + ylab("Proportion of Diseased Corals") + geom_smooth(se = F, col = "gray")
  
healthtemp <- ggplot(temp.dis, aes(x = tempC,  y = prop, group= Coral_stat)) + geom_point(size = 3, aes(color = Coral_stat)) + theme_bw() + xlab(templab) + ylab("Proportion of Corals") + geom_smooth(se = F, aes(color = Coral_stat), method = "lm")+ scale_color_manual(values = c("black","orange", "gray"), name = "Health", labels = c("Dead","Healthy", "Disease"))



plot_grid(templot, healthtemp, labels = "auto")
## `geom_smooth()` using formula 'y ~ x'

temp.prop.lm <- lm(prop~tempC*Coral_stat, data= temp.dis)
summary(temp.prop.lm)
## 
## Call:
## lm(formula = prop ~ tempC * Coral_stat, data = temp.dis)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.055072 -0.013532  0.002014  0.013404  0.053653 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -1.936543   0.266787  -7.259 1.57e-09 ***
## tempC              0.068640   0.008957   7.664 3.45e-10 ***
## Coral_statH        5.796832   0.377293  15.364  < 2e-16 ***
## Coral_statX        1.012797   0.377293   2.684  0.00963 ** 
## tempC:Coral_statH -0.179529   0.012667 -14.173  < 2e-16 ***
## tempC:Coral_statX -0.026392   0.012667  -2.084  0.04195 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02192 on 54 degrees of freedom
## Multiple R-squared:  0.988,  Adjusted R-squared:  0.9869 
## F-statistic:   891 on 5 and 54 DF,  p-value: < 2.2e-16
plot(temp.prop.lm)

car::Anova(temp.prop.lm, type = "III")
## Anova Table (Type III tests)
## 
## Response: prop
##                    Sum Sq Df F value    Pr(>F)    
## (Intercept)      0.025308  1   52.69 1.565e-09 ***
## tempC            0.028209  1   58.73 3.452e-10 ***
## Coral_stat       0.129381  2  134.68 < 2.2e-16 ***
## tempC:Coral_stat 0.112518  2  117.13 < 2.2e-16 ***
## Residuals        0.025937 54                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(temp.prop.lm)
##       (Intercept)             tempC       Coral_statH       Coral_statX 
##       -1.93654313        0.06864041        5.79683246        1.01279693 
## tempC:Coral_statH tempC:Coral_statX 
##       -0.17952941       -0.02639181
temp.prop.lmX <- lm(prop~tempC, data= temp.dis[which(temp.dis$Coral_stat == "X"),])
temp.prop.lmD <- lm(prop~tempC, data= temp.dis[which(temp.dis$Coral_stat == "D"),])
temp.prop.lmH <- lm(prop~tempC, data= temp.dis[which(temp.dis$Coral_stat == "H"),])

summary(temp.prop.lmX)
## 
## Call:
## lm(formula = prop ~ tempC, data = temp.dis[which(temp.dis$Coral_stat == 
##     "X"), ])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.055072 -0.013532  0.004643  0.018411  0.037225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.92375    0.30080  -3.071 0.006583 ** 
## tempC        0.04225    0.01010   4.184 0.000558 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02471 on 18 degrees of freedom
## Multiple R-squared:  0.493,  Adjusted R-squared:  0.4648 
## F-statistic:  17.5 on 1 and 18 DF,  p-value: 0.0005581
summary(temp.prop.lmD)
## 
## Call:
## lm(formula = prop ~ tempC, data = temp.dis[which(temp.dis$Coral_stat == 
##     "D"), ])
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0177993 -0.0058061  0.0007323  0.0074983  0.0166745 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.936543   0.130610  -14.83 1.57e-11 ***
## tempC        0.068640   0.004385   15.65 6.30e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01073 on 18 degrees of freedom
## Multiple R-squared:  0.9316, Adjusted R-squared:  0.9278 
## F-statistic:   245 on 1 and 18 DF,  p-value: 6.301e-12
summary(temp.prop.lmH)
## 
## Call:
## lm(formula = prop ~ tempC, data = temp.dis[which(temp.dis$Coral_stat == 
##     "H"), ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.04191 -0.01733  0.00134  0.01598  0.05365 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.86029    0.32556   11.86 6.12e-10 ***
## tempC       -0.11089    0.01093  -10.14 7.15e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02674 on 18 degrees of freedom
## Multiple R-squared:  0.8512, Adjusted R-squared:  0.8429 
## F-statistic: 102.9 on 1 and 18 DF,  p-value: 7.148e-09